Original Idea was to target JAMA - thinking one of the Research Letter article type:
Research Letters are concise, focused reports of original research. These should not exceed 600 words of text and 6 references and may include up to 2 tables or figures. Online supplementary material is only allowed for brief additional and absolutely necessary methods but not for any additional results or discussion. Research Letters may have no more than 7 authors. The text should include the full name, academic degrees, and a single institutional affiliation for each author and the email address for the corresponding author. Other persons who have contributed to the study may be indicated in an Acknowledgment, with their permission, including their academic degrees, affiliation, contribution to the study, and an indication if compensation was received for their role. Letters must not duplicate other material published or submitted for publication. In general, Research Letters should be divided into the following sections: Introduction, Methods, Results, and Discussion. They should not include an abstract, but otherwise should follow all of the guidelines in Manuscript Preparation and Submission Requirements. Letters not meeting these specifications are generally not considered.
In 2019, following the most devastating wildfire season in the state’s history, California’s largest utiility, Pacific Gas and Electric (PG&E) began using the practice of pre-emptive deenergization of utility line to avoid setting fires. During the 2019 Fire Season in California these events, refered to as Public Safety Power Shutoffs (PSPS), were used on several occasions.
When PSPS events occur residents experience a loss of power that can last beyond the
Spatial data of PSPS affected areas initially received from PG&E was cross-referenced with reports on their website.
There are some reports on PG&E’s website for which we do not have spatial data about affected areas.
Spatial data from PG&E contained boundaries for 12 separate areas, which - using layer names - were cross-referenced to reports for 7 specific de-energization events (shown on the map below, click for PG&E’s estimates start/end times and number of affected customers).
event1 <- filter(pge_events, event == "June_7-9")
event2 <- filter(pge_events, event == "September_25-27")
event3 <- filter(pge_events, event == "October_5")
event4 <- filter(pge_events, event == "October_9-12")
event5 <- filter(pge_events, event == "October_23-25")
event6 <- filter(pge_events, event == "October_26-November_1")
event7 <- filter(pge_events, event == "November_20-21")
pal2 <- colorQuantile(palette = c("#2172B4","#9DCAE1","#A4DBA1","#008836"), n=4, domain = 0:100)
maphpi <- merge(tracts_sp, {
hpi[,.(ct10, hpi2_pctile)]
})
leaflet() %>%
# addProviderTiles(providers$Esri.WorldStreetMap, group = "Street Map") %>%
# addProviderTiles(providers$Stamen.Terrain, group = "Terrain") %>%
# addProviderTiles(providers$Esri.WorldImagery, group = "Sattelite") %>%
addProviderTiles(providers$Stamen.Toner, group = "B/W") %>%
addPolygons(data = maphpi,
color = "#444444",
weight = 1,
smoothFactor = 0.1,
fillOpacity = 0.5,
fillColor = ~ pal2(maphpi$hpi2_pctile),
stroke = FALSE,
popup = paste0("This location is in the top ",round(maphpi$hpi2_pctile,0),"the percentile of all tracts in the state for the healthy places index"),
group = "HPI"
) %>%
addPolygons(
data = st_zm(event1),
color = "#252525",
weight = 2,
smoothFactor = 0.1,
fillOpacity = 0.7,
fillColor = "#377eb8",
stroke = TRUE,
label = ~event1$event,
popup = paste0("This PSPS event spanned ", round(event1$Duration_hrs,0), " hours \n (beginning ", as.character(event1$start_time), " and ending at ", event1$end_time,"). \n Ultimately PG&E estimates that it affected ",format(event1$customers, big.mark = ',')," customers."),
group = "June 7-9 (blue)"
) %>%
addPolygons(
data = st_zm(event2),
color = "#252525",
weight = 2,
smoothFactor = 0.1,
fillOpacity = 0.7,
fillColor = "#4daf4a",
stroke = TRUE,
label = ~event2$event,
popup = paste0("This PSPS event spanned ", round(event2$Duration_hrs,0), " hours \n (beginning ", as.character(event2$start_time), " and ending at ", event2$end_time,"). \n Ultimately PG&E estimates that it affected ",format(event2$customers, big.mark = ',')," customers."),
group = "Sep 25-27 (green)"
) %>%
addPolygons(
data = st_zm(event3),
color = "#252525",
weight = 2,
smoothFactor = 0.1,
fillOpacity = 0.7,
fillColor = "#984ea3",
stroke = TRUE,
label = ~event3$event,
popup = paste0("This PSPS event spanned ", round(event3$Duration_hrs,0), " hours \n (beginning ", as.character(event3$start_time), " and ending at ", event3$end_time,"). \n Ultimately PG&E estimates that it affected ",format(event3$customers, big.mark = ',')," customers."),
group = "Oct 5 (purple)"
) %>%
addPolygons(
data = st_zm(event4),
color = "#252525",
weight = 2,
smoothFactor = 0.1,
fillOpacity = 0.7,
fillColor = "#ffff33",
stroke = TRUE,
label = ~event4$event,
popup = paste0("This PSPS event spanned ", round(event4$Duration_hrs,0), " hours \n (beginning ", as.character(event4$start_time), " and ending at ", event4$end_time,"). \n Ultimately PG&E estimates that it affected ",format(event4$customers, big.mark = ',')," customers."),
group = "Oct 9-12 (yellow)"
) %>%
addPolygons(
data = st_zm(event5),
color = "#252525",
weight = 2,
smoothFactor = 0.1,
fillOpacity = 0.7,
fillColor = "#a65628",
stroke = TRUE,
label = ~event5$event,
popup = paste0("This PSPS event spanned ", round(event5$Duration_hrs,0), " hours \n (beginning ", as.character(event5$start_time), " and ending at ", event5$end_time,"). \n Ultimately PG&E estimates that it affected ",format(event5$customers, big.mark = ',')," customers."),
group = "Oct 23-25 (brown)"
) %>%
addPolygons(
data = st_zm(event6),
color = "#252525",
weight = 2,
smoothFactor = 0.1,
fillOpacity = 0.7,
fillColor = "#e41a1c",
stroke = TRUE,
label = ~event6$event,
popup = paste0("This PSPS event spanned ", round(event6$Duration_hrs,0), " hours \n (beginning ", as.character(event6$start_time), " and ending at ", event6$end_time,"). \n Ultimately PG&E estimates that it affected ",format(event6$customers, big.mark = ',')," customers."),
group = "Oct 26 - Nov 1 (red)"
) %>%
addPolygons(
data = st_zm(event7),
color = "#252525",
weight = 2,
smoothFactor = 0.1,
fillOpacity = 0.7,
fillColor = "#ff7f00",
stroke = TRUE,
label = ~event7$event,
popup = paste0("This PSPS event spanned ", round(event7$Duration_hrs,0), " hours \n (beginning ", as.character(event7$start_time), " and ending at ", event7$end_time,"). \n Ultimately PG&E estimates that it affected ",format(event7$customers, big.mark = ',')," customers."),
group = "Nov 20-21 (orange)"
) %>%
addLayersControl(
# baseGroups = c("Street Map", "Terrain", "Sattelite", "B/W"),
overlayGroups = c("HPI","June 7-9 (blue)", "Sep 25-27 (green)", "Oct 5 (purple)","Oct 9-12 (yellow)", "Oct 23-25 (brown)" , "Oct 26 - Nov 1 (red)", "Nov 20-21 (orange)"),
options = layersControlOptions(collapsed = FALSE)
)%>% # adds a button to return to the whole state view
addLegend("bottomleft",
pal = pal2,
values = na.omit(maphpi$hpi2_pctile),
labels = c("bottom 25%","","","healthiest 25%"),
opacity = 1,
title = "HPI Percentile") %>%
setView(lat = 39.302394, lng = -122.084003, zoom = 6) %>% #defaults the view of the original map to show the entire state
addEasyButton(easyButton(
icon="fa-globe", title="Zoom to State",
onClick=JS("function(btn, map){ map.setView([37.085206, -119.540085],5); }"))) %>%
hideGroup(c("June 7-9 (blue)", "Sep 25-27 (green)", "Oct 5 (purple)","Oct 9-12 (yellow)", "Oct 23-25 (brown)" , "Oct 26 - Nov 1 (red)", "Nov 20-21 (orange)")) # adds a button to return to the whole state view
# addLegend("bottomleft",
# pal = pal,
# values = factor(c("Moderate","High","VeryHigh"), levels = c("Moderate","High","VeryHigh")),
# opacity = 1,
# title = "Fire Hazard Severity Zones"
# ##lables = c("X%-X% (most vulnerable)","X%-X%","X%-X%","X%-X%","X%-X% (least"%. In this tract, the "mapTemp$def," is higher than "mapTemp),))##
# )In order to compare population affected and not affected by specific events, we query the American Community Survey 5-year estimates for several indicators of vulnerability.
% White non-Hispanic in the different areas by event.
############ DISABILITY ##############
race <- get_acs(
geography = "tract",
variables = c(
# numerator = "B01001A_001", # white alone
# numerator = "B01001B_001", # black alone
# numerator = "B01001C_001", # American Indian Alaska alone
# numerator = "B01001D_001", # Asian alone
# numerator = "B01001E_001", # API alone
# numerator = "B01001F_001", # other alone
# numerator = "B01001G_001", # 2 or more races
numerator = "B01001H_001", # white alone nonhispanic
# numerator = "B01001I_001", # Hispanic
denominator = "B01001_001" # total people
),
state = "CA"
) %>% data.table()
race <- race[, .(estimate = sum(estimate),
moe = moe_sum(estimate = estimate, moe = moe)),
by = .(GEOID, NAME, variable)]
try1 <- merge(race, {
readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
}, allow.cartesian = T)
# try1 <- na.omit(try1)
try1 <- try1[, se:= moe/1.645]
sums <- try1[,.(sum = sum(estimate, na.rm=T),
se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]
totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")
sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]
final <- merge(totals, sum_ses)
final <- final[, `:=` (proportion = numerator/denominator,
se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]
final <- final[, `:=` (LL = proportion - 1.96*se_prop,
UL = proportion + 1.96*se_prop)]
final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)
final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]
hc1 <- highchart() %>%
hc_chart(type = "bar") %>%
hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>%
hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
hc_add_theme(hc_theme_smpl()) %>%
hc_xAxis(title = list(text = "PSPS Event"),
type = 'category' ,
categories = c("June 7-9 (14hrs, 22k)", "Sep 25-27 (26hrs, 75k)","Oct 5 (18hrs, 11k) ","Oct 9-12 (90hrs, 729k)","Oct 23-25 (52hrs, 177k)","Oct 26-Nov 1 (129hrs, 941k)","Nov 20-21 (42hrs, 49k)")
) %>%
hc_yAxis(title = list(text = "% white (non-Hispanic) (with 95% CI)")) %>%
hc_title(text = "% white and non-Hispanic",
margin = 20, align = "left") %>%
hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
align = "left", style = list(fontWeight = "bold"))
hc1% Living in in the different areas by event.
############ DISABILITY ##############
poverty <- get_acs(
geography = "tract",
variables = c(
numerator = "B05010_002", # below 1.00 ratio
numerator = "B05010_010", # income ration 1.0 - 1.99
denominator = "B05010_001" # total workers
),
state = "CA"
) %>% data.table()
poverty <- poverty[, .(estimate = sum(estimate),
moe = moe_sum(estimate = estimate, moe = moe)),
by = .(GEOID, NAME, variable)]
try1 <- merge(poverty, {
readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
}, allow.cartesian = T)
# try1 <- na.omit(try1)
try1 <- try1[, se:= moe/1.645]
sums <- try1[,.(sum = sum(estimate, na.rm=T),
se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]
totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")
sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]
final <- merge(totals, sum_ses)
final <- final[, `:=` (proportion = numerator/denominator,
se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]
final <- final[, `:=` (LL = proportion - 1.96*se_prop,
UL = proportion + 1.96*se_prop)]
final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)
final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]
hc1 <- highchart() %>%
hc_chart(type = "bar") %>%
hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>%
hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
hc_add_theme(hc_theme_smpl()) %>%
hc_xAxis(title = list(text = "PSPS Event"),
type = 'category' ,
categories = c("June 7-9 (14hrs, 22k)", "Sep 25-27 (26hrs, 75k)","Oct 5 (18hrs, 11k) ","Oct 9-12 (90hrs, 729k)","Oct 23-25 (52hrs, 177k)","Oct 26-Nov 1 (129hrs, 941k)","Nov 20-21 (42hrs, 49k)")
) %>%
hc_yAxis(title = list(text = "% below 200% poverty level (with 95% CI)")) %>%
hc_title(text = "% below 200% poverty level",
margin = 20, align = "left") %>%
hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
align = "left", style = list(fontWeight = "bold"))
hc1% White non-Hispanic in the different areas by event.
############ DISABILITY ##############
insurance <- get_acs(
geography = "tract",
variables = c(
numerator = "B27001_005", # male under 6 no health insurance
numerator = "B27001_008", # male 6-18 no health insurance
numerator = "B27001_011", # male 19-25 no health insurance
numerator = "B27001_014", # male 26-34 no health insurance
numerator = "B27001_017", # male 35-44 no health insurance
numerator = "B27001_020", # male 45-54 no health insurance
numerator = "B27001_023", # male 55-64 no health insurance
numerator = "B27001_026", # male 65-74 no health insurance
numerator = "B27001_029", # male 75+ no health insurance
numerator = "B27001_033", # male under 6 no health insurance
numerator = "B27001_036", # male 6-18 no health insurance
numerator = "B27001_039", # male 19-25 no health insurance
numerator = "B27001_042", # male 26-34 no health insurance
numerator = "B27001_045", # male 35-44 no health insurance
numerator = "B27001_048", # male 45-54 no health insurance
numerator = "B27001_051", # male 55-64 no health insurance
numerator = "B27001_054", # male 65-74 no health insurance
numerator = "B27001_057", # male 75+ no health insurance
denominator = "B27001_001" # total people
),
state = "CA"
) %>% data.table()
insurance <- insurance[, .(estimate = sum(estimate),
moe = moe_sum(estimate = estimate, moe = moe)),
by = .(GEOID, NAME, variable)]
try1 <- merge(insurance, {
readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
}, allow.cartesian = T)
# try1 <- na.omit(try1)
try1 <- try1[, se:= moe/1.645]
sums <- try1[,.(sum = sum(estimate, na.rm=T),
se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]
totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")
sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]
final <- merge(totals, sum_ses)
final <- final[, `:=` (proportion = numerator/denominator,
se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]
final <- final[, `:=` (LL = proportion - 1.96*se_prop,
UL = proportion + 1.96*se_prop)]
final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)
final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]
hc1 <- highchart() %>%
hc_chart(type = "bar") %>%
hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>%
hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
hc_add_theme(hc_theme_smpl()) %>%
hc_xAxis(title = list(text = "PSPS Event"),
type = 'category' ,
categories = c("June 7-9 (14hrs, 22k)", "Sep 25-27 (26hrs, 75k)","Oct 5 (18hrs, 11k) ","Oct 9-12 (90hrs, 729k)","Oct 23-25 (52hrs, 177k)","Oct 26-Nov 1 (129hrs, 941k)","Nov 20-21 (42hrs, 49k)")
) %>%
hc_yAxis(title = list(text = "% without Health Insurance (with 95% CI)")) %>%
hc_title(text = "% without Health Insurance",
margin = 20, align = "left") %>%
hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
align = "left", style = list(fontWeight = "bold"))
hc1Below is an example for the % of people over 65 in the different areas by event.
old <- tidycensus::get_acs(geography = "tract",
variables = c(numerator = "B01001_020",
numerator = "B01001_021",
numerator = "B01001_022",
numerator = "B01001_023",
numerator = "B01001_024",
numerator = "B01001_025",
numerator = "B01001_044",
numerator = "B01001_045",
numerator = "B01001_046",
numerator = "B01001_047",
numerator = "B01001_048",
numerator = "B01001_049",
denominator = "B01001_001"),
state = "CA") %>% data.table()
old <- old[,.(estimate = sum(estimate),
moe = tidycensus::moe_sum(estimate = estimate, moe = moe)),
by= .(GEOID, NAME, variable)]
try1 <- merge(old, {
readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
}, allow.cartesian = T)
# try1 <- na.omit(try1)
try1 <- try1[, se:= moe/1.645]
sums <- try1[,.(sum = sum(estimate, na.rm=T),
se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]
totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")
sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]
final <- merge(totals, sum_ses)
final <- final[, `:=` (proportion = numerator/denominator,
se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]
final <- final[, `:=` (LL = proportion - 1.96*se_prop,
UL = proportion + 1.96*se_prop)]
final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)
final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]
hc1 <- highchart() %>%
hc_chart(type = "bar") %>%
hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>%
hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
hc_add_theme(hc_theme_smpl()) %>%
hc_xAxis(title = list(text = "PSPS Event"),
type = 'category' ,
categories = c("June 7-9 (14hrs, 22k)", "Sep 25-27 (26hrs, 75k)","Oct 5 (18hrs, 11k) ","Oct 9-12 (90hrs, 729k)","Oct 23-25 (52hrs, 177k)","Oct 26-Nov 1 (129hrs, 941k)","Nov 20-21 (42hrs, 49k)")
) %>%
hc_yAxis(title = list(text = "% over 65 (with 95% CI)")) %>%
hc_title(text = "% of population over 65",
margin = 20, align = "left") %>%
hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
align = "left", style = list(fontWeight = "bold"))
hc1% of people under age 5 in the different areas by event.
young <- tidycensus::get_acs(geography = "tract",
variables = c(numerator = "B01001_003", # male under 5
numerator = "B01001_027", # female under 5
denominator = "B01001_001"),
state = "CA") %>% data.table()
young <- young[,.(estimate = sum(estimate),
moe = tidycensus::moe_sum(estimate = estimate, moe = moe)),
by= .(GEOID, NAME, variable)]
try1 <- merge(young, {
readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
}, allow.cartesian = T)
# try1 <- na.omit(try1)
try1 <- try1[, se:= moe/1.645]
sums <- try1[,.(sum = sum(estimate, na.rm=T),
se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]
totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")
sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]
final <- merge(totals, sum_ses)
final <- final[, `:=` (proportion = numerator/denominator,
se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]
final <- final[, `:=` (LL = proportion - 1.96*se_prop,
UL = proportion + 1.96*se_prop)]
final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)
final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]
hc1 <- highchart() %>%
hc_chart(type = "bar") %>%
hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>%
hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
hc_add_theme(hc_theme_smpl()) %>%
hc_xAxis(title = list(text = "PSPS Event"),
type = 'category' ,
categories = c("June 7-9 (14hrs, 22k)", "Sep 25-27 (26hrs, 75k)","Oct 5 (18hrs, 11k) ","Oct 9-12 (90hrs, 729k)","Oct 23-25 (52hrs, 177k)","Oct 26-Nov 1 (129hrs, 941k)","Nov 20-21 (42hrs, 49k)")
) %>%
hc_yAxis(title = list(text = "% under 5 (with 95% CI)")) %>%
hc_title(text = "% of population under 5",
margin = 20, align = "left") %>%
hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
align = "left", style = list(fontWeight = "bold"))
hc1% of people with a disability in the different areas by event.
############ DISABILITY ##############
disability <- tidycensus::get_acs(geography = "tract",
variables = c(numerator = "C18130_010", # with a disability 18-64
numerator = "C18130_017", # with a disability 65+
denominator = "C18130_001"),
state = "CA") %>% data.table()
disability <- disability[,.(estimate = sum(estimate),
moe = tidycensus::moe_sum(estimate = estimate, moe = moe)),
by= .(GEOID, NAME, variable)]
try1 <- merge(disability, {
readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
}, allow.cartesian = T)
# try1 <- na.omit(try1)
try1 <- try1[, se:= moe/1.645]
sums <- try1[,.(sum = sum(estimate, na.rm=T),
se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]
totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")
sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]
final <- merge(totals, sum_ses)
final <- final[, `:=` (proportion = numerator/denominator,
se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]
final <- final[, `:=` (LL = proportion - 1.96*se_prop,
UL = proportion + 1.96*se_prop)]
final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)
final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]
hc1 <- highchart() %>%
hc_chart(type = "bar") %>%
hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>%
hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
hc_add_theme(hc_theme_smpl()) %>%
hc_xAxis(title = list(text = "PSPS Event"),
type = 'category' ,
categories = c("June 7-9 (14hrs, 22k)", "Sep 25-27 (26hrs, 75k)","Oct 5 (18hrs, 11k) ","Oct 9-12 (90hrs, 729k)","Oct 23-25 (52hrs, 177k)","Oct 26-Nov 1 (129hrs, 941k)","Nov 20-21 (42hrs, 49k)")
) %>%
hc_yAxis(title = list(text = "% with a disability (with 95% CI)")) %>%
hc_title(text = "% of population living with a disability",
margin = 20, align = "left") %>%
hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
align = "left", style = list(fontWeight = "bold"))
hc1% of Households receiving Food Stamps/SNAP in the past 12 months in the different areas by event.
############ DISABILITY ##############
snap <- get_acs(
geography = "tract",
variables = c(
numerator = "B22001_002", # receiving SNAP
# numerator = "B22001_003", # receiving SNAP over 60
denominator = "B22001_001" # total
),
state = "CA"
) %>% data.table()
snap <- snap[, .(estimate = sum(estimate),
moe = moe_sum(estimate = estimate, moe = moe)),
by = .(GEOID, NAME, variable)]
try1 <- merge(snap, {
readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
}, allow.cartesian = T)
# try1 <- na.omit(try1)
try1 <- try1[, se:= moe/1.645]
sums <- try1[,.(sum = sum(estimate, na.rm=T),
se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]
totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")
sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]
final <- merge(totals, sum_ses)
final <- final[, `:=` (proportion = numerator/denominator,
se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]
final <- final[, `:=` (LL = proportion - 1.96*se_prop,
UL = proportion + 1.96*se_prop)]
final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)
final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]
hc1 <- highchart() %>%
hc_chart(type = "bar") %>%
hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>%
hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
hc_add_theme(hc_theme_smpl()) %>%
hc_xAxis(title = list(text = "PSPS Event"),
type = 'category' ,
categories = c("June 7-9 (14hrs, 22k)", "Sep 25-27 (26hrs, 75k)","Oct 5 (18hrs, 11k) ","Oct 9-12 (90hrs, 729k)","Oct 23-25 (52hrs, 177k)","Oct 26-Nov 1 (129hrs, 941k)","Nov 20-21 (42hrs, 49k)")
) %>%
hc_yAxis(title = list(text = "% receiving food stamps (with 95% CI)")) %>%
hc_title(text = "% of Households receiving Food Stamps (past 12 months)",
margin = 20, align = "left") %>%
hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
align = "left", style = list(fontWeight = "bold"))
hc1% without vehicle available in the different areas by event.
############ DISABILITY ##############
vehicles <- get_acs(
geography = "tract",
variables = c(
# numerator = "B22001_002", # receiving SNAP
numerator = "B08141_002", # workers without a vehicle
denominator = "B08141_001" # total workers
),
state = "CA"
) %>% data.table()
vehicles <- vehicles[, .(estimate = sum(estimate),
moe = moe_sum(estimate = estimate, moe = moe)),
by = .(GEOID, NAME, variable)]
try1 <- merge(vehicles, {
readRDS("//mnt/projects/ohe/PSPS/PSPS_yn_master_tract.RDS")
}, allow.cartesian = T)
# try1 <- na.omit(try1)
try1 <- try1[, se:= moe/1.645]
sums <- try1[,.(sum = sum(estimate, na.rm=T),
se_sum = (sum(se^2, na.rm=T)^0.5)), by = .(variable, PSPS_yn, event)]
totals <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "sum")
sum_ses <- sums %>% dcast.data.table(formula = PSPS_yn + event ~ variable, value.var = "se_sum") %>% .[,.(PSPS_yn, event, SE_num = numerator, SE_pop = denominator)]
final <- merge(totals, sum_ses)
final <- final[, `:=` (proportion = numerator/denominator,
se_prop = (1/denominator)*(SE_num^2 - (((numerator/denominator)^2)*SE_pop^2))^0.5)]
final <- final[, `:=` (LL = proportion - 1.96*se_prop,
UL = proportion + 1.96*se_prop)]
final <- final[, event := factor(event, levels = c("June_7-9", "September_25-27","October_5","October_9-12","October_23-25","October_26-November_1","November_20-21"))] %>% setorder(event)
final <- final[, PSPS_yn := ifelse(PSPS_yn == "No","Not Affected by PSPS", "Affected by PSPS")]
hc1 <- highchart() %>%
hc_chart(type = "bar") %>%
hc_add_series(final, "point", hcaes(y = 100*proportion, x = event, group = PSPS_yn))%>%
hc_add_series(final, "errorbar", hcaes(low = 100*LL, high = 100*UL, x = event, group = PSPS_yn)) %>%
hc_add_theme(hc_theme_smpl()) %>%
hc_xAxis(title = list(text = "PSPS Event"),
type = 'category' ,
categories = c("June 7-9 (14hrs, 22k)", "Sep 25-27 (26hrs, 75k)","Oct 5 (18hrs, 11k) ","Oct 9-12 (90hrs, 729k)","Oct 23-25 (52hrs, 177k)","Oct 26-Nov 1 (129hrs, 941k)","Nov 20-21 (42hrs, 49k)")
) %>%
hc_yAxis(title = list(text = "% without a vehicle (with 95% CI)")) %>%
hc_title(text = "% of without access to a vehicle",
margin = 20, align = "left") %>%
hc_subtitle(text = "Tracts Impacted by PSPS vs Not Impacted",
align = "left", style = list(fontWeight = "bold"))
hc1